The instantaneous phase, amplitude, and frequency of a sawtooth wave

can changing the smoothness of a signal change the instantaneous freq dynamics? (i.e. more sharp transitions)



In [1]:

    
%config InlineBackend.figure_format = 'retina'
%matplotlib inline

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

import seaborn as sns
sns.set_style('white')

1. Simulate a 20Hz sawtooth wave



In [2]:

    
# Define sawtooth shape in some number of samples
x1 = np.array([0,.05,.2,1,.9,.8,.7,.6,.5,.4,.3,.2,.1,.05,.01])
t1 = np.arange(len(x1))

# Interpolate sawtooth so it has 50 samples (50ms = 20Hz wave)
from scipy import interpolate
f = interpolate.interp1d(t1, x1)
t2 = np.linspace(0,len(t1)-1,50)
x2 = f(t2)

# Tile the new sawtooth to last 5 seconds
x = np.tile(x2,100)
x = x - np.mean(x)
Fs = 1000.
t = np.arange(0,5,.001)

# Plot sawtooth
plt.figure(figsize=(5,2))
plt.plot(t, x)
plt.ylim((-.7,.7))
plt.xlim((0,.5))
plt.xlabel('Time (s)')
plt.ylabel('Voltage (a.u.)')









    Out[2]:





<matplotlib.text.Text at 0xab77470>

2. Filter in 13-30Hz band

Filter is very short (like steph said in email)



In [3]:

    
from misshapen import nonshape
x_filt, _ = nonshape.bandpass_default(x, (13,30),Fs, w=1.5,rmv_edge=False)

3. Calculate instantaneous measures



In [4]:

    
x_amp = np.abs(sp.signal.hilbert(x_filt))
x_phase = np.angle(sp.signal.hilbert(x_filt))



In [5]:

    
# Instantaneous freq
x_freq = np.diff(x_phase)
x_freq[x_freq<0] = x_freq[x_freq<0]+2*np.pi
x_freq = x_freq*Fs/(2*np.pi)

4. Visualize



In [6]:

    
samp_plt = range(1000,1200)

plt.figure(figsize=(5,5))
plt.subplot(3,1,1)
plt.plot(t[samp_plt],x[samp_plt],'k')
plt.plot(t[samp_plt],x_filt[samp_plt],'r')
plt.plot(t[samp_plt],x_amp[samp_plt],'b')
plt.xlim((t[samp_plt][0],t[samp_plt][-1]))
plt.xticks([])
plt.ylabel('raw (black)\nfiltered (red)\nInst. Amp. (blue)')
plt.subplot(3,1,2)
plt.plot(t[samp_plt],x_phase[samp_plt],'k')
plt.xlim((t[samp_plt][0],t[samp_plt][-1]))
plt.xticks([])
plt.ylabel('Inst. Phase (rad)')
plt.subplot(3,1,3)
plt.plot(t[samp_plt],x_freq[samp_plt],'k')
plt.xlim((t[samp_plt][0],t[samp_plt][-1]))
plt.xlabel('Time (s)')
plt.ylabel('Inst. Freq. (Hz)')









    Out[6]:





<matplotlib.text.Text at 0xc59f4e0>

5. Example with real data



In [7]:

    
x2 = np.load('C:/gh/bv/misshapen/exampledata.npy')



In [8]:

    
from misshapen import nonshape
x2_filt, _ = nonshape.bandpass_default(x2, (13,30),Fs, w=3,rmv_edge=False)

x_amp = np.abs(sp.signal.hilbert(x2_filt))
x_phase = np.angle(sp.signal.hilbert(x2_filt))

# Instantaneous freq
x_freq = np.diff(x_phase)
x_freq[x_freq<0] = x_freq[x_freq<0]+2*np.pi
x_freq = x_freq*Fs/(2*np.pi)



In [9]:

    
samp_plt = range(3000,5000)

plt.figure(figsize=(10,6))
plt.subplot(4,1,1)
plt.plot(t[samp_plt],x2[samp_plt],'k')
plt.plot(t[samp_plt],x2_filt[samp_plt],'r')
plt.plot(t[samp_plt],x_amp[samp_plt],'b')
plt.xlim((t[samp_plt][0],t[samp_plt][-1]))
plt.xticks([])
plt.ylabel('raw (black)\nfiltered (red)\nInst. Amp. (blue)')
plt.subplot(4,1,2)
plt.plot(t[samp_plt],x_phase[samp_plt],'k')
plt.xlim((t[samp_plt][0],t[samp_plt][-1]))
plt.xticks([])
plt.ylabel('Inst. Phase (rad)')
plt.subplot(4,1,3)
plt.plot(t[samp_plt],x_freq[samp_plt],'k')
plt.xlim((t[samp_plt][0],t[samp_plt][-1]))
plt.xticks([])
plt.ylim((0,30))
plt.ylabel('Inst. Freq. (Hz)')
plt.subplot(4,1,4)
plt.plot(t[samp_plt],x_amp[samp_plt],'k')
plt.ylabel('Inst. Amp.')
plt.xlabel('Time (s)')
ax = plt.gca()
ax2 = ax.twinx()
ax2.plot(t[samp_plt],x_freq[samp_plt],'r')
plt.xlim((t[samp_plt][0],t[samp_plt][-1]))
plt.ylim((0,30))
plt.ylabel('Inst. Freq. (Hz)',color='r')









    Out[9]:





<matplotlib.text.Text at 0xeabde80>

6. Plot relationship between inst. amp. and inst. freq.



In [15]:

    
plt.figure(figsize=(6,6))
plt.plot(x_freq, x_amp[1:],'k.',alpha=.01)
plt.xlim((0,100))
plt.xlabel('Inst. Freq. (Hz)')
plt.ylabel('Inst. Amp. (uV)')









    Out[15]:





<matplotlib.text.Text at 0x1180ad30>